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ABSTRACT 

We have investigated the influence of a velocity shear surface on the linear and non- 
linear development of the CD kink instability of force-free helical magnetic equilibria 
in 3D. In this study we follow the temporal development within a periodic computa- 
tional box and concentrate on flows that are sub-Alfvenic on the cylindrical jet's axis. 
Displacement of the initial force-free helical magnetic field leads to the growth of CD 
kink instability. We find that helically distorted density structure propagates along 
the jet with speed and flow structure dependent on the radius of the velocity shear 
surface relative to the characteristic radius of the helically twisted force-free magnetic 
field. At small velocity shear surface radius the plasma flows through the kink with 
minimal kink propagation speed. The kink propagation speed increases as the velocity 
shear radius increases and the kink becomes more embedded in the plasma flow. A de- 
creasing magnetic pitch profile and faster flow enhance the influence of velocity shear. 
Simulations show continuous transverse growth in the nonlinear phase of the instability. 
The growth rate of the CD kink instability and the nonlinear behavior also depend on 
the velocity shear surface radius and flow speed, and the magnetic pitch radial profile. 
Larger velocity shear radius leads to slower linear growth, makes a later transition to 
the nonlinear stage, and with larger maximum amplitude than occur for a static plasma 
column. However, when the velocity shear radius is much greater than the characteristic 
radius of the helical magnetic field, linear and non-linear development can be similar to 
the development of a static plasma column. 

Subject headings: instabilities - MHD - methods: numerical - galaxies: jets 



1. Introduction 



Relativistic jets occur in black hole binary star systems (microquasars) (e.g., Mirabel & 
Rodon'guez 1999), occur in active galactic nuclei (AGN) (e.g., Urry & Padovani 1995; Ferrari 1998; 
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Meier et al. 2001), can be associated with neutron stars and pulsar wind nebulae, e.g., Crab nebula 
jet (Weisskopf et al. 2000), and arc thought responsible for the gamma-ray bursts (GRBs) (e.g., 
Zhang & Mcszaros 2004; Piran 2005; Mcszaros 2006). It is thought that these jets arc powered and 
collimated hydromagnetically (e.g., Blandford 2000). Such Poynting flux dominated outflows arise 
on magnetic field lines threading the horizon of a rotating black hole (Blandford & Znajek 1977; 
Blandford & Payne 1982), and form as magnetic flux is pumped by the inflowing gas or generated 
by dynamo action in a surrounding accretion disk (e.g., Kudoh & Kaburaki 1996). The observed 
high degree of jet collimation is generally attributed to hoop stress from the toroidal magnetic field 
acting in concert with an external pressure. 

Jet generation simulations using general relativistic magnetohydrodynamics (GRMHD) (e.g., 
De ViUiers et al. 2003, 2005; Hawley & Krolik 2006, McKinney & Gammie 2004; McKinncy 2006; 
McKinney k Blandford 2009; Beckwith et al. 2008; Hardee et al. 2007; Komissarov &: Barkov 
2009; Penna et al. 2010) codes show development of magneto-rotational instability (MRI) (Balbus 
Sz Hawley 1998) and angular momentum transfer in the accretion disk, leading to diffusion of matter 
and magnetic field inwards, and unsteady outfiows near a centrifugally supported funnel wall. In 
general, GRMHD simulations with spinning black holes indicate jet production consisting of a 
Poynting-flux high Lorcntz factor spine with t; ~ c, and a matter dominated sheath with t; ~ c/2 
possibly embedded in a lower speed, v « c, disk/coronal wind. Circumstantial evidence such 
as the requirement for large Lorentz factors suggested by the TeV BL Lacs when contrasted with 
much slower observed motions (Ghisellini et al. 2005) suggests such a spine-sheath morphology, 
although alternative interpretations arc also possible (Georganopoulos & Kazanas 2003; Stern & 
Poutanen 2008; Bromberg &; Levinson 2009; Giannios et al. 2009). 

Strongly magnetized relativistic outflows are typically produced from rotating bodies (neutron 
stars, black holes and accretion disks). A toroidal magnetic field (i?^) is wound up in such outflows 
and in the far zone becomes dominant because the poloidal field (Bp) falls off faster with expansion 
and distance. In configurations with strong toroidal magnetic field, the current-driven (CD) kink 
mode is unstable. This instability excites large-scale helical motions that can strongly distort or even 
disrupt the system. For static cylindrical force-free equilibria, the well-known Kruskal-Shafranov 
criterion indicates that the instability develops if the length of the plasma column, I, is long enough 
for the field lines to go around the cylinder at least once (e.g., Bateman 1978): \Bp/B(f,\ < IjliiR. 
For relativistic force-free static configurations, the linear instability criteria have been studied by 
several authors (Istomin & Pariev 1994, 1996; Begelman 1998; Lyubarskii 1999; Tomimatsu et al. 
2001; Narayan et al. 2009). In a more realistic case, rotation and shear could significantly affect 
the instability criterion. 

Twisted structures are observed in many AGN jets on sub-parsec, parsec and kiloparsec scales 
(e.g., Gomez et al. 2001; Lobanov & Zensus 2001). Non-relativistic and relativistic simulations of 

magnetized jet formation and/or propagation have showed helical structures attributed to CD kink 
instability (e.g., Lery et al. 2000; Ouyed et al. 2003; Nakamura & Meier 2004; Nakamura et al. 
2007; MoU et al. 2008; Moh 2009; McKinney h Blandford 2009; Carey & Sovinec 2009). In the 
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absence of CD kink instability and resistive relaxation, helical structures may be attributed to the 
Kelvin-Helmholtz (KH) helical instability driven by velocity shear at the boundary between the jet 
and the surrounding medium (e.g., Hardee 2004, 2007) or triggered by precession of the jet ejection 
axis (Begelman et al. 1980). It is still not clear whether current driven, velocity shear driven or 
jet precession is responsible for the observed structures, or whether these different processes are 
responsible for the observed twisted structures at different spatial scales. 

In a previous paper we performed 3D relativistic MHD simulations and investigated the de- 
velopment of CD kink instability in a static plasma column with a force-free helical magnetic field 
(Mizuno et al. 2009a). We found that the initial configuration was strongly distorted but not 
disrupted by the kink instability. Although static configurations (or more generally rigidly mov- 
ing flows considered in the proper reference frame) are the simplest ones for studying the basic 
properties of the kink instability, in a realistic case, rotation and shear motions could significantly 
affect the CD kink instability. In this paper we investigate the influence of a velocity shear surface 
on the stability and nonlinear behavior of relativistic sub-Alfvenic flow. Because sub-Alfvenic flow 
is stable to the Kelvin-Helmholtz (KH) instability, we can focus on the development of CD kink 
instability. In this paper we describe the numerical method and setup used for our simulations in 
§2, present our results in §3, in §4 compare our results to instability expectations and conclude, 
and in the Appendix present multi-dimensional numerical code tests. 



2. Numerical Method and Setup 

In order to study time evolution of the CD kink instability in the relativistic MHD (RMHD) 
regime, we use the 3D GRMHD code "RAISHIN" in Cartesian coordinates. RAISHIN is based 
on a 3 + 1 formalism of the general relativistic conservation laws of particle number and energy- 
momentum. Maxwell's equations, and Ohm's law with no electrical resistance (ideal MHD con- 
dition) in a curved spacetime (Mizuno et al. 2006). In the RAISHIN code, a conservative, high- 
resolution shock-capturing scheme is employed. The numerical fluxes are calculated using the 
HLL approximate Riemann solver, and flux-interpolated constrained transport (flux-CT) is used 
to maintain a divergence-free magnetic fleld . The RAISHIN code performs special relativistic 
calculations in Minkowski spacetime by choosing the appropriate metric. The RAISHIN code has 
proven to be accurate to second order and has passed a number of numerical tests including highly 
relativistic cases and highly magnetized cases in both special and general relativity (Mizuno et al. 
2006). The results of multidimensional numerical test problems are shown in the Appendix. The 



^Constrained transport schemes are used to maintain divergence- free magnetic fields in the RAISHIN code. This 
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are interpolated to the cell center and as a result the scheme becomes non-conservative even though we solve the 
conservation laws (Komissarov 1999). 
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multidimensional numerical test problems show that the MC slope-limiter scheme performs best 
on most of the tests and we have used this scheme for reconstruction. 

In our simulations we choose a force-free helical magnetic field for the initial configuration 

(Mizuno ct al. 2009a). A force- free configuration is a reasonable choice for the strong magnetic 
field cases that we study here. In general, the force-free equilibrium of a static cylinder is described 
by the equation 

In particular we choose the following form for poloidal {B^) and toroidal (-B^) components of 
magnetic field determined in the laboratory frame 

R = (o\ 
' [1 + {R/afY ' 



Bo l [l + {R/aff--l-2a{R/af 

^ {R/a)[l + {R/aYY\ 2a - 1 ' ^ ' 

where R is the radial position in cylindrical coordinates normalized by a simulation scale unit 
L = \, Bq parameterizes the magnetic field amplitude, a is the characteristic radius of the column, 
and a is the pitch profile parameter. The pitch parameter, P = RB^/B^, determines the radial 
profile of the magnetic pitch, and larger indicates increasing helical pitch of the magnetic field 
lines. With our choice for the force-free field, the pitch parameter can be written as 



P - {R/a?]J ^ (i?/a)2]2a _ 1 _ 2a{R/af ' 

If the pitch profile parameter 0.5 < a < 1, the magnetic helicity increases with radius. When a = 1, 
the magnetic helicity is constant and if a > 1, the magnetic helicity decreases. This definition for 
the pitch parameter and the force-free helical field has been chosen to be the same as that used in 
previous non-relativistic work (Appl et al. 2000; Baty 2005) and in our previous relativistic work 
(Mizuno et al. 2009a) for purposes of comparison. For our modestly relativistic sub-Alfvenic flow 
speeds, the jet density and magnetic fields within the jet when determined in the jet flow frame 
are reduced slightly relative to their values in the laboratory frame by the flow Lorentz factor (e.g.. 
Anile 1989; Komissarov 1999). 

The simulation grid is periodic along the axial z direction. As a consequence the allowed axial 
wavelengths are restricted to A = Lz/n < L^, with n a positive integer and is the grid length. 
The grid is a Cartesian {x,y,z) box of size 4L x 4L x 3L with grid resolution of AL = L/40. 
The grid resolution is the same in all directions. In simulations we choose a characteristic radius, 
a = (1/4)L. In terms of a, the simulation box size is 16a x 16a x 12a and the allowed axial 
wavelengths are restricted to A = 12a/n < 12a. We impose outflow boundary conditions on the 
transverse boundaries at x = y = ±2L (zb8a). This simulation grid is the same as that used for 
case B in Mizuno et al. (2009a). We checked the influence of grid resolution for the case of a 
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static plasma column using four different grid resolutions from 20 to 60 computational zones per 
simulation length unit L = 8a in Mizuno et al. (2009a). We found that growth does depend on grid 
resolution, and our choice of AL = L/40 is sufficient to capture the growth of CD kink instability. 

We consider a low gas pressure medium with constant p = Po = 0.02 in units of poc^ for the 
equilibrium state, and a non-uniform density profile decreasing with the magnetic field strength as, 
p = piB^ with pi = 6.25/90- We choose a density decreasing oc in order to keep the Alfven speed 
inside the velocity shear surface above the flow speed. The equation of state is that of an ideal gas 
with p = (T — l)pe, where e is the specific internal energy density and the adiabatic index F = 5/3 
[§. The specific enthalpy is /i = 1 + e/c^ -\-p/ pc? . The magnetic field amplitude is Bq = 0.4 in units 
of a/ 47rpoc^ leading to a low plasma-/3 near the axis. The sound speed is c^/c = {Tp/ph)^/"^ and 
the Alfven speed is given by va/c = [B'^/{ph + B'^)]^/'^. 

We choose two different jet velocities: (case s) slow, Vj = 0.2c and (case f) fast, Vj = 0.3c. 
In order to study the effect of the velocity, we perform the simulations with four different velocity 
shear surface radii: Rj = 1/8L (a/2) , 1/4L (a), 1/2L (2a), IL (4a). Results are compared to 
those for a static plasma column (no flow) as the reference. We also investigate the effect of 
different magnetic pitch profiles with: (case CP) constant pitch, a = 1, and (case DP) decreasing 
pitch, a = 2.0. The radial profiles of the magnetic field components, the magnetic pitch, and the 
sound and Alfven speeds determined in the lab frame for the different cases are shown in Figure 
1. For all cases the sound and Alfven speeds on the axis are c^q = 0.178 c and vao = 0.364 c, 
respectively, when determined in the lab frame. We have chosen constant and decreasing pitch 
cases only because previous simulation studies of the CD kink instability of a static plasma column 
show that increasing magnetic pitch is not too different from the constant pitch case (Mizuno et al. 
2009a). The jet velocity is sub-Alfvenic inside the velocity shear radius for almost all cases with the 
exception of the constant pitch case with the largest shear radius Rj = 4a and fastest jet velocity 
Vj = 0.3c. In this paper we focus on the development of CD kink instability as the sub-Alfvenic jet 
is stable to the velocity shear driven Kelvin-Helmholtz (KH) instability. 

To break the symmetry the initial MHD equilibrium configuration is perturbed by a radial 
velocity in all region with profile given by 



The amplitude of the perturbation is 5v = 0.01c with radial width Rp = 0.5L (2a), and we choose 
m = 1 and n = 1. This is identical to imposing {m,n) = (—1,-1), because of the symmetry 



^ This adiabatic index assumes the plasma is cold, Cs ^ c. This condition is generally fulfilled in our simulations 
even though the sound speed is greater than the relativistic limit of c/y/S at large radial position (see Fig. If). 
The region where the CD kink instability occurs is near the axis and in this region an adiabatic index F = 5/3 is 
appropriate. We have checked the dependence of our simulation results on different adiabatic indicies (see Appendix 
A in Mizuno et al. 2009a). We found no significant difference in results for different adiabatic indicies. 
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Fig. 1. — Radial profile of (a) the toroidal magnetic field B^, (b) the axial magnetic field i?^, (c) the 
Alfven speed, va/c, (d) magnetic and gas (dash-dot line) pressure, (e) the rest mass density, p, and 
(f) the sound speed, Cg/c, all determined in the lab frame. Solid lines indicate the constant pitch 
case and dashed lines indicate the decreasing pitch case. Dotted lines in (c) show the locations of 
the velocity shear radius, Rj/a = 0.5, 1.0, 2.0, and 4.0 at Vj = 0.2c. 



between (m, n) and (— m, —n) pairs. The various different cases that we have considered are listed 
in Table 1. 



3. Results 

3.1. Constant Helical Pitch: Vj = 0.2 c 

Figure 2 shows the time evolution of a density isosurface for constant helical pitch with Vj = 
0.2c and Rj = 2a (CPs2a) where the time, t, is in units of tc = L/c = Aa/c (light travel time 
across the largest velocity shear surface radius, Rj = 4a, considered in this study). Displacement of 
the initial force-free helical magnetic field by growth of the CD kink instability leads to a helically 
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Case: CPs2a (constant pitch, Vj^O.lc, Rj^2a) 




Fig. 2. — Time evolution of three-dimensional density isosurfaces with a transverse slice at z = 
for case CPs2a. The time, t, is in units of = L/c. Color shows the logarithm of the density with 
solid magnetic field lines. Velocity vectors are shown by the arrows. 

twisted magnetic filament wound around the density isosurface. In the nonlinear phase, helically 
distorted density structure shows continuous transverse growth and propagates in the flow direction. 
This propagation of the helical kink structure does not occur for a static plasma column (Mizuno 
et al. 2009a). 



In order to investigate the dependence of kink growth and propagation on the location of the 
velocity shear surface, we consider four velocity shear radii from Rj = a/2 to 4a (see Table 1). The 
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effect of different radii on the growth of the CD kink for constant magnetic pitch and Vj = 0.2c is 
shown in Figure 3. In Figure 3 we show the time evolution of the volume-averaged kinetic, Ei^in,xyj 

Case: CPs (constant pitch, Vj=0.2c) 




0. 0070 1- -I 

10 20 30 40 50 60 70 
time 

Fig. 3. — Time evolution of (a) Ekin,xy and (b) Emag,xy for constant pitch (a = 1.0) with Vj = 0.2c 
at different jet radii: Rj/a = 0.5 (red solid line), 1.0 (orange dashed line), 2.0 (green dotted line) 
and 4.0 (blue dash-dotted line). For reference a static plasma column case (no jet) is also shown 
as black dash double dotted lines. 

and magnetic, Emag,xy, energy transverse to the z-axis determined within a cylinder of radius 
R/L < 1.0 {R < 4a) as an indicator of the growth of the CD kink instability. The wavelength 
of the kink is A = 12a. Note that according to the Kruskal-Shafranov criterion, the instability 
develops at A > 27ra. The instability growth rate reaches a maximum at Xmax ~ 10a, the exact 
coefficient being dependent on the transverse distribution of the density and magnetic pitch, and 
also possibly the location and magnitude of the velocity shear. For the case of constant pitch 
and uniform density Appl et al. (2000) found Xmax = 8.43a and a corresponding growth rate of 
^max = 0.133v AO /cL- In general, one can use the estimate T^ax ~ O.Ivao/cl in the rest frame of the 
kink. For a moving kink we might expect the temporal growth rate in the lab frame to be reduced 
with r oc where is the moving kink Lorentz factor, e.g., Narayan et al. (2009). Change in 
the evolution of Eki^^xy and Emag,xy indicate an initial linear growth phase at t < (35 — 55)fc with 
duration depending on the velocity shear radius, and followed by a nonlinear evolution phase. In 
all cases, the initial growth phase is characterized by an exponential increase in Ei^in,xy by about 3 
orders of magnitude to a maximum amplitude followed by a slow decline in the nonlinear phase. By 
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fitting the linear portion of the slope in Ekin,xy between the amplitudes of 10~^ and 5 x 10~^ we can 
determine an e-folding time where Tg = At/ In 50 and At is the time interval. The e-folding times 
can be found in Table 2 in §4. In general, the e-folding time increases as the velocity shear radius 
increases for a/2 < Rj < 2a. The e-folding time at Rj = 4a decreases but is still significantly longer 
than for the static plasma column. The time evolution trend of Emag,xy is opposite to the time 
evolution of Ekin,xy Emag,xy gradually decreases in the early linear growth phase, then exhibits an 
initial rapid decrease into the nonlinear phase to a minimum followed by a slight increase at later 
times. 

At the smallest velocity shear radius the behavior of Ekin,xy and Eynag,xy are very similar to 
that of a static plasma column. Increased difference in behavior from that of a static plasma column 
appears as the velocity shear radius increases to Rj = a and 2a. For these CctSGS ctS the velocity 
shear radius increases the growth rate of the CD kink slows, and Ekin,xy achieves a somewhat higher 
maximum amplitude with a later transition to the nonlinear stage. As the shear radius increases 
Emag,xy exhibits a more gradual decline in the transition to the nonlinear stage than for the static 
plasma column. However, for the largest velocity shear radius, Rj = 4a, the difference relative to 
the static plasma column in the linear and early non-linear phase is reduced, and the growth rate of 
the CD kink instability is faster and with lower maximum E^in^xy than for the cases with velocity 
shear radius Rj = a and Rj = 2a. But note that Ekin,xy becomes slightly larger than the static 
case at the longest comparable simulation times. 

When the velocity shear radius is much larger than the characteristic radius of the helical 
magnetic field, we would expect the growth of CD kink instabihty to approach that of a static 
plasma column moving with respect to the observer. However, in the flow reference frame the 
Alfven speed for the sub-Alfvenic jet is not the same as that for a static plasma column in addition 
to a small relativistic clock effect. Thus, we do not expect the initial growth rate of largest jet 
radius case to perfectly match that of the static plasma column as determined in the observer 
(simulation) rest frame. The effect of the velocity shear is greatest for the case with velocity shear 
radius Rj = 2a but also significant when Rj = a and 4a. 

Figure 4 shows a density isosurface at i = 50 for the constant pitch cases with Vj = 0.2c 
at different velocity shear radii. In the linear growth phase, the behavior of the growing kink is 
almost the same for different velocity shear radius. However, in the nonlinear phase, the behavior 
of the kink is different for different velocity shear radius. For the smallest radius, Rj = a/2, the 
kink does not propagate along the z-axis significantly. Transverse amplitude growth dominates this 
case and is very similar to that of a static plasma column. In general the flow appears to follow 
the helical twist indicated by the density isosurface. Flow velocities in the transverse x-y plane 
are larger than flow along the z-axis. When Rj = a flow in the x-y plane is reduced somewhat 
relative to flow along the z-axis and the kink propagates slowly along the z-axis. As the radius 
increases the kink propagates along the z-axis more rapidly (see more detail in Figure 5) while 
continuing transverse growth. Flow is less twisted helically as the velocity shear radius increases 
and at Rj = 4a we see a helical kink embedded within and moving with the flow. Recall that the 
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time evolution of the volume averaged transverse kinetic energy, i.e., Fig. 3, were the most similar 
to a static plasma column for smallest and largest shear radius. The principal difference between 
the smallest and largest velocity shear radius cases lies in the flow morphology relative to the kink 
morphology indicated by the density isosurface. 



Case: CPs (constant pitch, Vj=0.2c) 
Rj=l/2a Rj=a 




Fig. 4. — Three-dimensional density isosurface at t = 50 with a transverse slice at z = for constant 
helical pitch with Vj = 0.2c at different velocity shear radii: (a) Rj/a = 0.5, (b) 1.0, (c) 2.0, and 
(d) 4.0. Color shows the logarithm of the density with velocity vectors indicated by the arrows. 
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In Figure 5 we follow the time evolution of the x and y— position of the maximum density 
in the xy plane at z = 1.5L = 6a. Transverse growth of the kink is revealed by displacement of 
this maximum away from zero in the xy plane. In all cases significant displacement of the density 
maximum begins at i > 20, when growth is still within the linear regime. An oscillation of the 
maximum indicates rotation around the z-axis in the xy plane and is related to the propagation 
speed of the kink. Since the kink wavelength is A = 12a in all simulations, the kink propagation 
speed is = Xv = [12a/(At x 4a)]c = (3/At)c where At is the length of time for one rotation in 
simulation time units. 

Case: CPs (constant pitch, Vj=0.2c) 
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Fig. 5. — Time evolution of (a) x and (b) y-position of maximum density in the xy plane at 
z/L = 1.5 (6a) for constant pitch (a = 1.0) with Vj = 0.2c at different shear radii: Rj/a = 0.5 (red 
solid line), 1.0 (orange dashed line), 2.0 (green dotted line) and 4.0 (blue dash-dotted line). The 
static plasma column is shown as a black dash-double dotted line. 

At the smallest shear radius (red solid lines), no oscillation is evident and transverse growth 
is very similar to that of the static plasma column (black dash double dotted lines). The lack of 
measurable propagation implies that the flow moves through the growing helical twist at nearly the 
flow speed, as suggested by the velocity vectors in panel (a) of Figure 4. When the shear radius 
is Rj = a (orange dashed lines), only a partial oscillation occurs by the end of the simulation. 
Our best estimate of At/ 2 > 25 comes from panel (b) in Figure 5 with a maximum and minimum 
displacement of pm,y occurring at t ~ 44 and t > 69, respectively. This implies a kink propagation 
speed of ffc/c < 0.06. The kink is propagating along jet axis while transverse growth continues but 
the propagation speed is slow. In this case the flow is still helically twisted but with less than the 
apparent helical twist of the density isosurface shown in panel (b) of Figure 4. 
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More than one complete oscillation is evident for a shear radius Rj = 2a (see the green dotted 
lines) and visual inspection suggests that the oscillation period becomes longer with time. At t < 45 
three measurements indicate {At/2) ~ 10 ± 1, where "±" indicates the range of the measurements. 
This implies a kink propagation speed Wfe/c ~ 0.15 ± 0.01. At t > 40 three measurements indicate 
{At/ 2) ~ 15 ±3. This implies a kink propagation speed w^/c ~ 0.10 it 0.02. The kink is propagating 
rapidly along the axis but as transverse growth continues the kink slows. In this case the flow is 
only modestly helically twisted and with much less than the apparent helical twist of the density 
isosurface shown in panel (c) of Figure 4. 

Multiple oscillations are evident for a shear radius Rj = 4a (see blue dash-dotted lines) and 
again the oscillation period becomes longer with time. Five measurements at t < 45 provide 
{At/2) ~ 7.5 ± 0.5, and five measurements at t > 45 provide {At/2) --^ 8.5 ± 1. Here the initial 
kink speed ffc/c ~ 0.20 ± 0.01 is approximately equal to the flow speed. At later times the kink 
speed tifc/c ~ 0.16 ± 0.02 has slowed significantly. At early time the smaller transverse amplitude 
of the kink is embedded within the flow and the flow exhibits little helical twist. The oscillation 
amplitude, indicative of the transverse amplitude of the kink, increases up to about t ~ 45 and 
then declines slightly with time. This suggests a decrease in transverse amplitude growth as the 
growing kink encounters the velocity shear surface and slows. 

We see from these results that kink propagation and flow morphology are strongly dependent 
on the velocity shear radius relative to the characteristic radius of the force-free plasma column. 
At one extreme the fluid flows through the kink with helicity comparable to that of the kink and at 
the other extreme the kink is embedded within a more uniform flow. The kink propagation speeds 
are listed in Table 2 in §4. 

3.2. Decreasing Helical Pitch: Vj = 0.2 c 

In Figure 6 we show the time evolution of the volume-averaged kinetic, E]^in,xyi ^■iid magnetic, 
Emag,xy, energy transverse to the ^;-axis determined within a cylinder of radius R/L < 1.0 {R < 4a), 
e.g.. Figure 3. Here we see an initial linear growth phase, t < (25 — 30)tc, and subsequent nonlinear 
evolution phase. The exponential increase in Ekin^xy hy about 3 orders of magnitude to a maximum 
amplitude followed by a slow decline in the nonlinear phase and behavior of Emag,xy is similar 
to what was found for the constant pitch cases. By fitting the linear portion of the slope in 
Ekin,xy between the amplitudes of 10~® and 5 x 10""^ we can determine an e-folding time where 
Te = At/ In 50 and At is the time interval. The c- folding times can be found in Tabic 2 in §4. In 
this case, the e-folding time is increased relative to the static plasma column by about the same 
amount only for velocity shear radii 2a < Rj < 4a. However, the decreasing helical pitch results 
in more rapid growth in the linear phase, and makes a transition to the nonlinear phase in about 
60% of the time for the constant pitch cases shown in Figure 3. This more rapid growth is similar 
to results for a static plasma column (Mizuno et al. 2009a). Thus, the growth rate trends found 
for the static plasma column are maintained in the presence of sub-Alfvenic flow. The maximum 
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Case: DPs (decrease pitch, Vj=0.2c) 




0. 0035 r j 

10 20 30 40 50 60 70 
time 

Fig. 6. — Time evolution of (a) Ekin,xy and (b) Emag,xy for decreasing helical pitch (a = 2.0) with 
Vj = 0.2c at different shear radii: Rj/a = 0.5 (red solid line), 1.0 (orange dashed line), 2.0 (green 
dotted line) and 4.0 (blue dash-dotted line). The static plasma column case (no flow) is also shown 
as black dash-double dotted lines. 

amplitude of Ekin,xy is comparable to that found for constant helical pitch. However, the amplitude 
of Ef:in,xy declines less in the non-linear phase relative to the static case than was found for the 
constant pitch cases. We note however that the constant pitch cases might have exhibited similar 
behavior at longer timescales. Here the decline in Ekin,xy is clearly the least when Rj = 2a with 
decline increasing for Rj = Aa, a and a/2, respectively. The same trend is evident in Figure 3 for 
constant helical pitch. 

In this set of simulations we follow the development of the kink to a much longer time relative 
to the transition time from the linear growth phase to the non-linear phase. Again we find that 
the largest effects occur for velocity shear radius Rj = 2a but now confirm relatively large effects 
for Rj = a and Rj = 4a. However, even for Rj = a/2 we clearly see the infiuence of velocity shear 
when compared with the results for constant pitch. We conclude that non-linear development is 
more infiuenced by a velocity shear surface in the case of decreasing helical pitch. 

Figure 7 shows a density isosurface for the decreasing pitch cases with Vj = 0.2c for different 
velocity shear radius at t = 50. The decreasing pitch cases shown here all appear similar to results 
shown for static plasma columns (Mizuno et al. 2009a). In the linear growth phase the properties 
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Case: DPs (decrease pitch, Vj=0.2c) 
Rj=l/2a Rj=a 




Fig. 7. — Three-dimensional density isosurface at t = 50 with transverse shces at z = for decrease 
pitch with Vj = 0.2c of different jet radius (a) Rj/a = 0.5, (b) 1.0, (c) 2.0, and (d) 4.0. Color shows 
the logarithm of the density with velocity vectors indicated by the arrows. 

are almost same for different shear radius but after transition to the nonlinear phase the behavior 
of the kink is different for different shear radius. In general, the behavior is similar to what was 
found for constant helical pitch. For the smallest radius, Rj = a/2, the kink does not propagate 
along the z-axis significantly. In general the flow appears to follow the helical twist indicated by 
the density isosurface. Flow velocities in the transverse x-y plane are larger than flow along the 
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z-axis. When Rj = a flow in the x-y plane is reduced somewhat relative to flow along the z-axis 
and the kink propagates slowly along the z-axis. As the radius increases the kink propagates along 
the z-axis more rapidly (see more detail in Figure 8) while continuing transverse growth. Flow is 
less twisted helically as the velocity shear radius increases. However, we note that when Rj = 4a 
we see considerably more indication of flow helicity than for the constant pitch case. 

In Figure 8 we follow the time evolution of the x and y— position of the maximum density in 
the xy plane at z = 1.5L = 6a, e.g.. Figure 5 for constant pitch cases. In all the decreasing pitch 
cases significant displacement of the density maximum begins at t < 15, when growth is still within 
the linear regime, and significant displacement occurs earlier than for the constant pitch cases. 



Case: DPs (decrease pitch, Vj=0.2c) 
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Fig. 8. — Time evolution of (a) x and (b) y-position of maximum density in the xy plane at 
z/L = 1.5 (6a) for decreasing pitch (a = 2.0) with Vj = 0.2c oat different shear radii: Rj/a = 0.5 
(red solid line), 1.0 (orange dashed line), 2.0 (green dotted line) and 4.0 (blue dash-dotted line). 
The static plasma column case (no flow) is shown as black dash-double dotted lines. 

At the smallest shear radius (red solid lines), no oscillation is evident and transverse growth 
is very similar to that of the static plasma column (black dash-double dotted lines). The lack of 
measurable propagation implies that the flow moves through the growing helical twist at nearly the 
flow speed, as suggested by the velocity vectors in panel (a) of Figure 7. When the shear radius 
is Rj = a (orange dashed lines), only a partial oscillation occurs by the end of the simulation. 
Our best estimate of At/4 ^ 20 comes from panel (a) in Figure 7 with a minimum and zero 
displacement of pm,x occurring at t ~ 30 and t ^ 50, respectively. This implies a kink propagation 
speed of Wfc/c ~ 0.04. Comparison to the comparable constant pitch case suggests a reduction in 
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the propagation speed but our level of accuracy does not make this a firm conclusion. The kink 

is propagating along jet axis while transverse growth continues but the propagation speed is slow. 
In this case the flow is still helically twisted but with less than the apparent helical twist of the 
density isosurface shown in panel (b) of Figure 7. 

More than one complete oscillation is evident for a shear radius Rj = 2a (sec the green dotted 
lines) and visual inspection suggests that the oscillation period becomes longer with time. At t < 35 
two measurements indicate (At/2) ~ 12.5±0.5, where "±" indicates the range of the measurements. 
This implies a kink propagation speed Vk/c ~ 0.12 it 0.005. At t > 35 two measurements indicate 
(At/2) ~ 24 ± 1. This implies a kink propagation speed Vk/c ~ 0.0625 ± 0.0025. The kink is 
propagating along the axis but as transverse growth continues the kink slows. Here the kink speed 
both at early and late times is significantly less than that found for the comparable constant pitch 
case. 

Multiple oscillations are evident for a shear radius Rj = 4a (sec blue dash-dotted lines) and 
again the oscillation period becomes longer with time. Four mcasTircments at t < 35 provide 
(At/2) ~ 7.5 ± 0.5, and three measurements at t > 35 provide (At/2) ~ 14.5 ± 1. Here the initial 
kink speed Vk/c ~ 0.20 ± 0.01 is approximately equal to the flow speed. At later times the kink 
speed Vk/c^ 0.105 ± 0.01 has slowed significantly. At early time the smaller transverse amplitude 
of the kink is embedded within the flow and the flow exhibits little helical twist. The oscillation 
amplitude, indicative of the transverse amplitude of the kink, increases up to about t ~ 50. Here 
we see transverse amplitude growth continuing to longer times and achieving larger amplitude than 
the comparable constant pitch case. The kink speed at later times is significantly less than that 
found for the comparable constant pitch case. 

We see from these results that kink propagation and flow morphology are again strongly 
dependent on the velocity shear radius relative to the characteristic radius of the force-free plasma 
column. In general, the presence of a velocity shear surface influences the propagation speed more 
than was found for the constant pitch cases. The slower kink propagation speeds found for the 
decreasing pitch cases are likely the result of the faster amplitude growth of the kink to a larger 
transverse amplitude. Thus, at the larger velocity shear radii the growing kink approaches the 
velocity shear surface more rapidly and/or more closely. The kink propagation speeds are listed in 
Table 2 in §4. 

3.3. Constant Helical Pitch: Vj = 0.3 c 

Figure 9 shows the time evolution of the volume-averaged kinetic energy {Ekin,xy) and magnetic 
energy {Ejnag,xy) transverse to the z-axis within a cylinder of radius R/L < 1.0 for constant pitch 
with Vj = 0.3c, similar to Figures 3 and 6. Change in the evolution of Ekin,xy and Efnag,xy indicate 
an initial linear growth phase at t < (35 — 60)tc with duration of the linear growth phase depending 
on the velocity shear radius. In all cases, the initial growth phase is characterized by an exponential 
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increase in Ei^in,xy by about 3 orders of magnitude to a maximum amplitude followed by a slow 
decline in the nonlinear phase. By fitting the linear portion of the slope in Ei^i^^xy between the 
amplitudes of 10~^ and 5 x 10~^ we can determine an e-folding time. The e-folding times are listed 
in Table 2 in §4. In general, the e-folding time increases as the velocity shear radius increases for 
a/2 < Rj < 2a. The e-folding time at Rj = 4a decreases but is still significantly longer than for 
the static plasma column. Evolution is qualitatively similar but quantitatively different from the 
constant pitch slower jet cases. In particular, the maximum value for Ekin,xy is about a factor of 
two larger than found for the constant pitch cases with vj = 0.2c, and is likely the result of the 
initial fiow kinetic energy being about a factor of two higher. As was found for the slower jet cases, 
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Fig. 9. — Time evolution of (a) Ekin,xy and (b) Emag,xy for constant pitch (a = 1.0) with Vj = 0.3c 
at different shear radii: Rj/a = 0.5 (red solid line), 1.0 (orange dashed line), 2.0 (green dotted 
line) and 4.0 (blue dash-dotted line). The static plasma column case (no flow) is shown as black 
dash-double dotted lines. 

the growth rate is slower and transition to the nonlinear phase occurs later as the shear radius 
increases from Rj = a/2 (red lines) to Rj = 2a (green dotted lines). As was found previously, the 
largest effects of velocity shear appear when Rj = 2a with smaller but still significant effects when 
Rj = a. 

The Rj = Aa simulation terminated before the maximum amplitude in E^i^^xy was achieved, 
but it is clear that the maximum amplitude in Ef^in^xy will occur at significantly later time than 
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when Rj = a. This result is different from the comparable slower flow case. This difference is likely 
the result of relativistic effects and will be considered further in §4. The three-dimensional helical 
structure of these faster jet cases is qualitatively similar to that of slower jet cases with constant 
pitch. Density isosurfaces, magnetic field lines and velocity vectors appear similar to those shown 
in Figures 2 &; 4. In Figure 10 we show the time evolution of the x and y— position of the maximum 
density in the xy plane at z = 1.5L = 6a, e.g., Figure 5 for the slower fiow constant pitch cases. 
Here significant displacement of the density maximum begins at t ~ 20, when growth is still within 
the linear regime. This is similar to the slower fiow constant pitch cases. At the smallest shear 



Case: CPf (constant pitch, Vj=0.3c) 




10 20 30 40 50 60 70 

time 



Fig. 10. — Time evolution of (a) x and (b) y-position of maximum density in the xy plane at 
z/L = 1.5 (6a) for constant pitch (a = 1.0) with vj = 0.3c at different jet radii: Rj/a = 0.5 (red 
solid line), 1.0 (orange dashed line), 2.0 (green dotted line) and 4.0(blue dash-dotted line). The 
static plasma column case (no fiow) is also shown as black dash-double dotted lines. 

radius (red solid lines), there is a suggestion of an oscillation in Pm,y (panel b) with a maximum 
negative displacement at t ~ 40 and At/2 > 30, although no oscillation is evident in pm,x (panel a). 
We can use this apparent oscillation to set an upper limit to the propagation speed of f^/c < 0.05. 
Thus, the fiow moves through the growing helical twist at nearly the fiow speed, as was found for 
the previous cases. 

When the shear radius is Rj = a (orange dashed lines), a full oscillation occurs and visual 
inspection suggests that the oscillation period increases with time. At t < 35 two measurements 
indicate {At/2) ~ 12.75 it 1, where "it" indicates the range. This implies a kink propagation speed 
Vk/c ~ 0.12 ± 0.01. At t > 35 two measurements indicate {At/2) ~ 24.75 it 2. This implies a kink 
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propagation speed v^/c ~ 0.06 ± 0.005. The kink is propagating along the axis but as transverse 
growth continues the kink slows. Here the kink speed at late times is similar to that found for the 
comparable slower flow speed constant pitch case. 

Multiple oscillations arc evident for a shear radius Rj = 2a (see the green dotted lines) and 
visual inspection suggests that the oscillation period becomes slightly longer with time. At t < 45 
measurements indicate {At/2) ~ 6.5 ± 0.5, where "it" indicates the range. This implies a kink 
propagation speed Vk/c ~ 0.23±0.015. At t > 45 measurements indicate {At/2) ~ 8.75±0.5. This 
implies a kink propagation speed Vk/c ~ 0.17 ± 0.01. The kink is propagating relatively rapidly 
along the axis at early times but as transverse growth continues the kink slows. 

Multiple oscillations are evident for a shear radius Rj = 4a (see blue dash-dotted lines) but 
there is no evidence that the oscillation period becomes longer with time for t < 53 when the 
simulation terminated. The multiple oscillations provide (At) ~ 9.8 it 0.4. Here the kink speed 
■Wfc/c ^ 0.305 ±0.01 is approximately equal to the flow speed. The oscillation amplitude, indicative 
of the transverse amplitude of the kink, increases up to the end of this simulation. This suggests 
that the kink remains embedded in the flow up to the end of this simulation. 

These results arc similar to those found for the slower flow cases with constant pitch. When 
the shear radius is larger, the propagation speed of the helical kink becomes faster. Quantitative 
comparison with individual slower jet cases (Fig. 5), on average indicates a more rapid propagation 
speed for each shear radius. The kink propagation speeds are listed in Table 2 in §4. 

4. Summary and Discussion 

We have investigated the development of the CD kink instability of a force-free helical magnetic 
field with a sub-Alfvenic velocity shear surface located at various radii relative to the characteristic 
radius of the magnetic field. We restricted this study to sub-Alfvenic shear as this regime is 
appropriate to the magnetically dominated flows thought to exist in the acceleration and collimation 
regions of relativistic jets. In this magnetically dominated parameter regime the flow is stable to 
the velocity shear driven Kelvin-Helmholtz instability so that we could focus solely on the effect 
of the shear flow on growth and propagation of the current driven kink and the velocity flow field 
accompanying the helically twisted kink. 

The growth of CD kink instability in the initial exponential growth phase is slower than found 

for a static plasma column. In general, the reduction in the growth rate is larger for constant 
magnetic pitch than for decreasing magnetic pitch, and is larger if the velocity is larger. In all 
cases, effects resulting from the presence of a velocity shear surface are largest when the velocity 
shear surface lies at or not too far outside, Rj = a &z, 2a, the characteristic radius, a, of the force-free 
magnetic field. For the slower speed, Vj = 0.2 c, and constant pitch case it was clear that the initial 
growth was least affected when the velocity shear surface was inside, Rj = a/2, or far outside, 
Rj = 4a, the characteristic radius. For the higher speed, Vj = 0.3 c, and constant pitch case initial 
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growth was again least affected when Rj = a/2 but relativistic effects slowed the initial growth 
when Rj = 4a. 

Transition to the nonlinear stage occurs at a later time but with larger maximum amplitude 
in the volume averaged transverse kinetic energy, Ekin,xy, as the velocity shear radius is increased 
for Rj < 2a, when compared to a static plasma column. However, when the velocity shear radius 
is far from the characteristic radius, Rj = 4a, the maximum amplitude is comparable to that of 
the static plasma column although that maximum is reached after longer time. In the absence of 
relativistic effects which slow the observed rate of growth, it is clear that the presence of a velocity 
shear surface has the strongest influence on both the linear and non-linear behavior of the growing 
kink when the velocity shear surface is located at or not too far outside the characteristic radius. 
In general, decreasing the magnetic pitch or increasing the flow velocity enhances the influence of 
the velocity shear surface. 

The location of the velocity shear surface has profound consequences for kink propagation 
and the associated flow field. For the velocity shear surface well inside the characteristic radius, 
transverse growth is similar to the static plasma column. In this case the plasma flows through 
the growing helical kink. For the velocity shear surface well outside the characteristic radius, the 
kink is embedded within and moves with the flow until the kink amplitude becomes large and the 
kink approaches the velocity shear surface. In this case the initial transverse growth is similar to 
that of a static plasma column advected with the flow and with growth computed in the proper 
reference frame. As the growing helically twisted kink approaches the velocity shear surface the 
kink slows and the flow field becomes slightly helically twisted. For velocity shear radii on the order 
of the characteristic radius there is a more intimate interaction between the growing kink and the 
flow field. In general, the kink propagates more slowly than the flow and slows as the amplitude 
increases. Thus, the flow field becomes more helically twisted as the kink amplitude increases. For 
these cases the flow helicity remains less than kink helicity. 

The Lorentz factors for slower and faster flows are 7 ~ 1.02 and 1.05, respectively. The 
growth rate of the CD kink instability depends on the Alfven velocity (e.g., Appl ct al. 2000) 
and also depends on relativistic time dilation. In the flow reference frame, the toroidal and axial 
magnetic field components are reduced by the Lorentz factor and the mass density is reduced by 
the Lorentz factor squared. The Alfven velocity is decreased slightly in the flow reference frame 
because h = 1 + e/c^ + P/pc? is slightly larger. Therefore in faster flow cases, we would expect the 
growth rate to be reduced primarily by time dilation related to the Lorentz factor of the moving 
kink, 7fc. Our results for kink e- folding times in the linear growth phase are summarized in Table 
2 along with the kink propagation speeds for the various different simulations. 

Quantitative comparison between our results and instability predictions is difficult because 
no sufficiently general stability analysis has been performed for magnetically dominated jets. The 
stability analysis performed by Narayan et al. (2009) considers the case of a "... rigid impenetrable 
wall at the outer cylindrical radius, -Rj.", that unfortunately is not appropriate to our simulations. 
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Nevertheless, we can make some comparison with previous results for static plasma columns and 
consider the implications for the spatial development of the instability. In all cases, the initial 
axisymmetric structure is strongly distorted by the kink instability, even though not disrupted. In 
general, the addition of a sub-Alfvenic velocity shear surface leads to slower temporal development 
of the instability than for the static case. Comparison between the static and moving kink constant 
pitch cases shows that the e-folding times (see Table 2) are significantly longer, up to ^ 16% 
and ~ 50% longer for slow and fast velocity shear, when the shear surface is located at twice the 
characteristic radius of the plasma column. When the velocity shear surface is located far outside 
the characteristic radius of the plasma column the kink is advected at about the flow speed in the 
linear growth regime and the e-folding times are shortened somewhat to about ~ 8% and ~ 40% 
longer than for the comparable static case. 

With the kink moving with the flow frame we would expect the growth rate in the flow frame 
to be related to the growth rate in the lab frame by the Lorentz factor of the flow. However, 
comparison between our static plasma column kink e-folding times and our moving kink e-folding 
times is complicated by length contraction in addition to time dilation. Time dilation increases 
the e-folding times from fluid to lab frames by the Lorentz factor. Length contraction means that 
our simulation box imposes a wavelength that appears longer in the fluid frame than in the lab 
frame by the Lorentz factor. In either frame the wavelength is longer than the fastest growing 
wavelength. Examination of our previous static case numerical results for the growth of E^i^^y 
at wavelengths of A = 12a and 16a (see Figure 2 in Mizuno et al. (2009a) for constant pitch 
cases A and B) indicates that the e-folding time increases about 7.5% faster than proportional 
to the wavelength over this wavelength range. Thus, length contraction means that our e-folding 
time observed for a static kink should convert approximately to our c-folding time observed for a 
moving kink by r™^ ~ 1.0757^r|*. The e-folding times for cases CPs4a and CPf4a when compared 
to the comparable static case, CPO, are significantly longer than the predicted 4% {vk ~ 0.2 c) 
and 11% {vk ~ 0.3 c) increase. We can only assume that a much larger velocity shear radius is 
required to further reduce the e-folding times for the propagating kinks to that predicted. On the 
other hand, this also means that temporal kink growth is significantly slowed even for a velocity 
shear surface at four times the characteristic radius. It is interesting to note that the fluid inertia 
increases by 7^ and if the growth time is also increased by the fluid inertia then the e-folding times 
in the lab frame should be increased by t^'" ^ 1.0757^r|* = 1.16 and 1.31 relative to the static 
case, and this increase comes much closer to the observed increase. We speculate that the large 
increase in e-folding times measured in the lab frame is partly a result of the increased inertia of 
the relativistically moving fluid. 

The characteristic time for the instability to affect strongly the initial structure varies from 
(25 — 30)tc for the decreasing pitch case to (35 — 60)tc for the constant pitch cases. For the constant 
pitch cases the characteristic time is roughly r ~ lOrg, with values for being dependent on 
the structure of the undisturbed state. In a jet context our perturbations remain static or can 
propagate with the flow frame depending on the location of the velocity shear surface. In order 
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to check whether the instabihty would affect a jet flow, one has to compare r with a propagation 
time. If wc identify Tg with the fastest growing wavelength, our present results suggest a scaling 
like T ~ 107^r|* with 3 > a > 1 for a moving kink and with velocity shear surface a few times 
the characteristic radius. Here a = 1 would correspond to time dilation only and a = 3 to time 
dilation plus inertial effects from the relativistically moving fluid. In this case the condition for the 
instability to affect the jet structure might be written as 

z > i^VkiA^) , (6) 

where we set lOrf = A{a/c) and < Vk < vj is a function of Rj/a and is sensitively dependent 
on the location of the velocity shear surface provided Rj/a « 10. This result suggests that the 
characteristic scale for kink development could be longer or very much shorter than for a kink 
simply advected with a broad flow for which z > 1007^ a. 

In order to find a more general criteria one has to know how the characteristic radius, a, and 
Lorentz factor, 7^, increase with distance. If the jet is narrow enough so that Qa^/c < z, where 
is the angular velocity at the base of the jet, one can use the scaling (Tchekhovskoy et al. 2008; 
Komissarov et al. 2009; Lyubarsky 2009) ^ Q,a/c and assume that 7^ = €7^. In this case one 
finds that the criterion for the kink instability can be written as 

^c>^(e^)"[l-(e^)-2]V'a. (7) 

The instability could affect the jet structure only if the jet expands slowly enough and/or the 
kink moves slowly enough. Assuming the parabolic shape for the jet, J7a/c = ^(Oz/c)*^ = ^x'^) 
where k < 1 and ^ ~ 1 are dimensionless numbers, one finds that the instability develops only if 
< l/(a + 1). In this case the characteristic scale for the development of the instability can be 
written as 

nz/c = X ~ [1 _ (,^^fc)-2]l/2[l-(a+l)fc] 

For e = 1, a = 1 and = Vj ^ c we recover the case for a kink advected with the flow field, eq.(9) 
in Mizuno et al. (2009a), ^Iz/c ~ (lO^^/^^"^''^ for A = 100. Here as 1 > e ^ l/^x^ Vk/vj « 1 
and ^z/c < 1 so the characteristic scale for development of instability can be very short. For a > 1, 
i.e., potential inertial effects, the characteristic time is lengthened and this has the potential for 
lengthening the characteristic scale for the development of the instability. This shows up in the 
constraint on k < l/(a + 1) < 1/2, i.e., only if the jet expands slowly enough. 

The 3D relativistic jet generation simulation performed by McKinney & Blandford (2009) 
indicates relatively rapid, less than 100 gravitational radii, but non-disruptive kink development 
over 500 gravitational radii. Our previous and present simulations for static and moving kinks 
suggest that the rapid but non-disruptive kink development in the jet generation simulation could 
be a result of a velocity shear surface located less or on order of the characteristic magnetic radius 
and a density increasing with radius. This combination would result in a slowly moving kink, hence 
rapid initial spatial development, but with non-linear growth slowed by the density increase and 
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accompanying Alfven speed decline with radius, increasing the Alfven crossing time and slowing 
spatial development. A non-linearly stabilizing increasing density profile is what might be expected 
for a Poynting flux jet core confined within a sensor slowly moving sheath, as appears to be predicted 
by jet generation simulations. Of course, a proper investigation of spatial growth requires stability 
simulations designed to study spatial kink development using more realistic flow, magnetic, and 
density profiles. 

In this paper, we considered sub-Alfvenic jet flow and focused on the development of the CD 
kink instability. If the velocity shear is super-Alfvenic, the flow can be KH unstable and CD 
unstable. Baty k, Keppens (2002) have investigated the interaction between KH and CD driven 
instabilities of a magnetized force-free cylindrical configuration via 3D MHD simulations in the 
non-relativistic regime. They found that the CD unstable modes provided a stabilizing effect on 
KH instability driven vortical structure. However, they assumed a relatively weak magnetic field 
in their simulations and studied the super-Alfvenic regime where the KH instability grows faster 
than the CD kink instability (e.g., Appl et al. 2000; Baty 2005). If the magnetic field is strong but 
the velocity shear is weakly super-Alfvenic, the growth rate of the KH instability can be less than 
or comparable to that of the CD instability. Even in the supcr-Alfvcnic regime the KH instability 
can be suppressed if a jet is embedded in a slower moving magnetized sheath that reduces the 
velocity shear to being effectively sub-Alfvenic (Hardee 2007; Mizuno et al. 2007). Such a sheath 
may exist around jets in the acceleration and collimation region. Thus, investigation of the weakly 
super-Alfvenic parameter regime will be particularly important to understanding the development 
of the twisted structures that are observed on relativistic jets. In future work we will investigate the 
spatial development of CD instability and the coupling between CD instability and KH instability. 
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A. Multidimensional Numerical Tests 

In this section we first summarize previous code results against known solutions and second 
we check our numerical simulation code results against those obtained by other codes using the 
same test problems. In RMHD there are not many problems with known solutions because a closed 
solution for the general RMHD Riemann problem has not yet been found. 

The code has been verified against a number of test and physical problems with known so- 
lutions. In Mizuno et al. (2006), a ID linear Alfven wave propagation test and a ID magnetized 
Bondi accretion flow test showed second-order convergence for our simulation code, and relativistic 
MHD shock-tube tests showed that our simulation code correctly handles the wave structure of 
both shocks and rarefactions even in an extreme relativistic regime. The code has been used to 
model several ID and 2D physical problems with known solutions. The code was successfully used 
to solve the ID Riemann problem for the deceleration of an arbitrarily magnetized relativistic flow 
injected into a static unmagnetized medium (Mizuno et al. 2009b). In particular, this physical 
problem involved comparing conditions for existence of the reverse shock against known relativistic 
shock jump conditions. The code has also been successfully used to solve in 2D the effect of mag- 
netic fields on an HD/MHD boost mechanism proposed by Aloy & Rezzolla (2006). The results to 
this physical problem involving a relativistic flow that flows parallel to an overpressure generated 
shock and rarefaction wave combination agreed with and extended the previous HD findings to 
MHD (Mizuno et al. 2008). 

We also have used our simulation code for several multidimensional (3D) physical problems 
involving comparison with known theoretical predictions. The code has been used to study spatial 
development of the Kelvin-Helmholtz instability on relativistic super- Alfvenic and trans- Alfvcnic 
cylindrical jets, and results (see Mizuno et al. 2007) were compared to predictions made by a linear 
stability analysis (Hardee 2007). The code has also been used to study temporal development 
of the current driven instability of a relativistic static plasma column, i.e., with magnetic energy 
density comparable to the plasma energy density including the rest mass energy, and results were 
successfully compared to linear stability analysis predictions and to previous non-relativistic nu- 
merical results (see Mizuno et al. 2009a). In both completely different instability regimes, the code 
delivered results that agreed with stability analysis predictions and with previous numerical results. 
These successful comparisons verify the code against known solutions to physical multidimensional 
problems. 

In general, we have not yet shown results for multidimensional test problems. In what follows 
we check our simulation code against multidimensional test problems using the MC slope-limiter 
scheme which we have used for the simulations in this paper. In all test simulations we use the HLL 
approximate Riemann solver to calculate numerical fluxes and the flux-CT scheme to maintain a 
divergence-free magnetic field. 
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A.l. Advection of a Magnetic Field Loop 

We consider the advection of a weak magnetic field loop in an initially uniform velocity field. 
In thermal pressure dominance the loop is transported as a passive scalar. The 2D test that we 
perform is a relativistic version of the magnetic loop advection problem. This tests the dissipative 
properties of the numerical scheme and the correct discretization balance of multidimensional terms 
(Gardiner &: Stone 2005, 2008; Mignone et al. 2010). A successful test results in the preservation 
of the initial loop. 

Following non-relativistic MHD tests (Gardiner & Stone 2005, 2008; Mignone et al. 2010), we 
employ a periodic computational box defined as —0.3 < x < 0.3 and —0.15 < y < 0.15 discretized 
on Nx X Nx/^ computational zones. Density and gas pressure are initially constant, and set p = 1.0 
and p = 0.328 respectively. The sound speed Cg = 0.3c for an adiabatic index V = 5/3. The 
magnetic field is defined through its magnetic vector potential as 

{ao + a2r2 if < r < 
Ao{R -r) if Ri<r<R, (Al) 
ii r>R, 

where Aq = 10"^, R = 0.09, Ri = 0.2R, 03 = -0.5Ao/Ri, oq = Ao{R - Ri) - a2Rj, and 
r = \J — y^. This modification of the vector potential in the r <R\ region, with respect to the 
original version of the test problem performed by Gardiner h Stone (2005), is done to remove the 
singularity in the loop's center that can cause spurious oscillations and erroneous evaluations of the 
magnetic energy (Mignone et al. 2010). We perform a case with no advection and an advection case 
in which the velocity of the flow is set up as Vx = 0.6c, Vy = 0.3c and Vz = 0.0c. The simulations 
are evolved until t = 1 when the loop has crossed through the periodic boundaries and returned to 
the center of the grid in the advection case. 

In Figure 11 the magnetic energy density (i?^) and magnetic file lines are shown for the no 
advection case (upper) and advection case (lower) with = 256. In the no advection case, the 
circular shape of the loop is perfectly preserved. There is no visible distortion or dissipation around 
the center and outer boundaries of the field loop. In the advection case, the circular shape of the 
loop is relatively well preserved although some flattening of the magnetic field lines is seen at left 
upper and right lower outer boundaries to the field loop. The magnetic energy density image shows 
internal structure comparable to that found in the Newtonian limit by Gardiner & Stone (2005) 
using first and second order CTU+CT schemes (Fig. 3 and 6 in Gardiner & Stone 2005). Mignone 
et al. (2010) performed a similar 2D magnetic loop advection test in the Newtonian limit using 
higher order schemes (third and fifth order) . Not surprisingly, the higher order schemes do a better 
job of maintaining uniformity in the magnetic energy density and in the loop structure. 

A quantitative measure of the magnetic field dissipation rate is indicated by the time evolution 
of the volume averaged magnetic energy density normalized to its initial value as shown in Figure 
12. Higher numerical resolution leads to a less diffusive result with less than 2% magnetic energy 



-26- 




Fig. 11. — Gray scale images of the magnetic energy density {B'^){left) and magnetic field lines 
(right) for the 2D magnetic field loop problem with no flow advection (upper) and with advection 
(lower) at t = 1 for Nx = 256. 
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Fig. 12. — Time evolution of the volume-averaged magnetic field energy density normalized to its 
initial value in the 2D field loop advection problem using the MC scheme with = 512 (solid), 
256 (dotted), and 128 (dashed). 




loss at the highest numerical resolution. Our result using = 128 finds almost the same magnetic 
energy dissipation, about 7% magnetic energy loss, as that found by Gardiner & Stone (2005) (see 
t = 1 in Fig. 7 of Gardiner &: Stone 2005) in the Newtonian limit using a comparable second- 
order scheme. Our quantitative results compare favorably with the higher order schemes tested by 
Mignone et al. (2010). Here our result using = 128 (~ 7% loss) is similar to their result using 
a third-order scheme (^^ 6.5% loss) at the same resolution (see i = 1 in Fig.A.9 of Mignone et al. 
2010). Our result using higher resolution = 256 shows less magnetic energy dissipation (~ 3% 
loss) at time t = 1. 



The three dimensional version of this problem is particularly challenging and see Gardiner & 

Stone (2008) for comparable second-order scheme results in the Newtonian limit. Correct evolution 
depends on how accurately the divergence- free condition is preserved and how well multidimensional 
MHD terms are balanced. The test consists of a computational box defined as —0.15 < x < 0.15, 
—0.15 <y< 0.15 and —0.3 < ^; < 0.3 discretized on x Nx/2 x A^^, computational zones with 

periodic boundary condition in all directions. Following the two-dimensional case, density and gas 
pressure arc initially constant, and we set p = 1.0 and p = 0.328. We show a no advection case 
and an advection case in which the velocity of the flow is set as {vx,Vy,Vz) = (0.3c, 0.3c, 0.6c). As 
for the two-dimensional case the vector potential is used to initialize the magnetic field in the 
coordinate system {xi,X2,X3) which is related to computational coordinate system {x,y,z) via the 
rotation 

xi = 00370; — sin 72; 
X2 = y 

X2 = sin ^x + cos 72;, 

where 7 = tan~^ 1/2. The initial vector potential is given by Eq. (Al), with r = \^x'f + x^- 

Figure 13 shows three-dimensional isovolumc images of the magnetic energy density for the no 
advection (left) and advection case (right) with A^^' = 256. In the no advection case, the cylindrical 
shape of field loops is perfectly preserved. Some small dissipation is seen around the center and 
outer boundaries of the field loops. In the advection case, the cylindrical shape of the field loops 
is relatively well preserved although there is some distortion and dissipation around the inner and 
outer boundaries of the field loops. This result is similar to that seen in the 2D field loop advection 
case (see Fig. 11) and appears comparable to the results shown in Figure 2 in Gardner k, Stone 
(2008). Mignone et al. (2010) performed a similar 3D magnetic loop advection test in the Newtonian 
limit using higher order schemes (third and fifth order). Not surprisingly, the higher order schemes 
do a better job preserving the loop structure and also display sharper boundaries. 

A quantitative measure of the magnetic field dissipation rate in 3D field loop advection problem 
is indicated by the time evolution of the volume averaged magnetic energy density normalized to 
its initial value as shown in Figure 14. As found for the 2D magnetic advection problem, higher 
numerical resolution leads to a less diffusive result with ^ 3% magnetic energy loss at the highest 
numerical resolution of Nx = 256. Our result using Nx = 128 finds a slightly greater magnetic 
energy dissipation (< 7% loss) when compared to the ~ 5% magnetic energy loss found by Gardiner 
&; Stone (2008) (see Fig. 1 in Gardiner &; Stone 2008) in the Newtonian limit using a comparable 
second-order scheme. Our magnetic energy loss using Nx = 128 (< 7% loss) is greater than that 
found by Mignone et al. (2010) (see Fig.A.9 in Mignone et al. 2010) using a third-order scheme 
(> 3% loss) at the same resolution. However, our result using higher resolution = 256 shows 
less magnetic energy dissipation (~ 3% loss). 

In summary, we obtain results very similar to previous non-relativistic results obtained using 
comparable second-order schemes (Gardiner & Stone 2005; Gardiner k. Stone 2008). At present 
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there are no published 2D or 3D relativistic loop test results for comparison. We conclude that 
our simulation code successfully passes both 2D and 3D test problems for the second-order MC 
slope-limiter and flux-CT schemes. We note that test simulations without the flux-CT scheme fail 
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Fig. 14. — Time evolution of the volume-averaged magnetic field energy density normalized to its 
initial value in the 3D field loop advection problem using the MC scheme with Nx = 256 (solid) 
and 128 (dashed). 

to pass these test problems. 

A. 2. 2D Cylindrical Explosion 

The cylindrical explosion test problem consists of a strong shock propagating into a magneti- 
cally dominated medium. Results for different cylindrical explosion problems in RMHD have been 
reported by several authors (e.g., Dubai 1991; van Putten 1995; Komissarov 1999; Del Zanna et 
al. 2003; Leismann et al. 2005; Mignone & Bodo 2006; Anton et al. 2010). Here, unlike the 2D 
loop test we can compare our results to identical relativistic test results obtained by other RMHD 
second-order codes. This test verifies that the generation of a high Lorentz factor flow is handled 
correctly. A successful test results in a dramatic difference in pressure and density between the 
ambient gas and the explosion zone. 

We have chosen a setup following that in Leismann et al. (2005) which is very similar to that 
in Komissarov (1999): a cylinder with high pressure (pc = 1), density (pc = 10"^), and radius 
0.8 is located in the center of a square Cartesian grid which initially contains a uniform, strong 
magnetic field. Between a radius of 0.8 to 1.0 the density and pressure smoothly decrease to those 
of a homogeneous ambient medium (p = 10~^ and p = 5 X 10~^). Initially, the magnetic field is 
in the x— direction, = 0.1, and the velocity is zero everywhere. The simulations are carried out 
until t = 4.0 on grid resolutions of 400^ and 800^, spanning a region of [—6, 6]^. 

Figure 15 shows the results from the 2D cylindrical explosion test problem att = 4.0. The outer 
fast shock has an almost circular shape. The innermost region is also almost circular and bounded 
by a reverse fast shock. Between these two fast shocks there are two density shells bounded on 
the outside by compressed magnetic field. Our test results are qualitatively identical to the results 
found by other second-order RMHD codes (see Figs. B3 & B4 in Leismann et al. 2005; see Fig. 4 
in Del Zanna et al. 2008). 
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Fig. 15. — 2D images of (a) density, (b) gas pressure, (c) Lorentz factor, and {d) magnetic pressure 
{Pm = B'^/2) of the 2D cylindrical explosion problem at t = 4.0 with 400^ resolution. 



In order to evaluate the quantitative difference between our simulation code and others, one- 
dimensional gas pressure and magnetic pressure profiles along the y axis are shown in Figure 16. 
The location of discontinuities is identical to the results in Leismann et al. (2005) and Del Zanna 
et al. (2007). The maximum value of the gas pressure is slightly larger than in Del Zanna et 
al. because of different numerical resolution but the same as in Leismann et al. at comparable 
numerical resolution. Slight differences in the one-dimensional magnetic pressure profile around 
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the maximum result from differences in numerical diffusivity in the different codes. 
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Fig. 16. — ID gas pressure (left) and magnetic pressure (right) profiles of the 2D cylindrical 
explosion problem at t = 4.0 along the y axis using the MC scheme with 400^ resolution. 

In summary, we obtain qualitatively identical and quantitatively almost identical results when 
compared to previous second order RMHD simulation code results. We conclude that our simulation 
code passes this test problem successfully. 

A.3. 3D Rotor 

Del Zanna et al. (2003) adapted the two-dimensional rotor problem from non-relativistic MHD 
performed by Balsara & Spicer (1999) and Toth (2000) to the relativistic case. Here we consider 
a three dimensional version of the relativistic rotor problem (Mignone et al. 2009). In a successful 
test the complicated pattern of shocks and torsional Alfven waves launched by the rotor is handled 
correctly. 

The initial condition consists of a sphere with radius rgp = 0.1 centered at the origin of the 
computational box taken to be the unit cube [—0.5,0.5]^. The sphere is heavier (psp = 10) than 
the surrounding medium (p = 1) and rapidly rotates around the z axis with an angular velocity 
LOsp = 9.95, i.e., with velocity components {vx,Vy,Vz) = C0sp{—y,x,0). The gas pressure, magnetic 
field and adiabatic index are constant everywhere, p = 1, B = (1,0,0), and T = 5/3. Exploiting 
the point symmetry, wc have carried out simulations until t = 0.4 at grid resolutions of 128^ and 
256'^. This test problem has resolution-dependent complexity since the maximum Lorentz factor in 
the initial set up depends on the grid resolution. 

Figure 16 shows two-dimensional density images in the xy plane at z = (perpendicular to the 
rotation axis) and the xz plane at y = (including the rotation axis) for the 256^ grid resolution 
case, dX t = 0.4. When the sphere starts rotating, torsional Alfven waves propagate outwards. 
The initial spherical structure collapses into a disk-like structure in the equatorial plane {z = 0) 
which generates shock waves propagating the vertical direction that can be seen in Figure 16. In 
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Fig. 17. — 2D density image of the 3D rotor problem at t = 0.4 (a) in the xy plane at z = and 
(&) in the xz plane at y = using the MC scheme with 256^ resolution. 

the xy plane, surrounding matter is pushed into a thin elliptical shell enclosed by a tangential 
discontinuity. In Mignone et al. (2009) the thin shell shows a distinct octagonal-like shape when 
they use a five wave HLLD Riemann solver scheme to calculate numerical fluxes that makes a 
less diffusive transition at rotational (Alfven) discontinuities. The whole structure is embedded 
in a radially expanding spherical fast rarefaction wave front. In the xz plane our spindle shaped 
structure appears nearly identical to the results shown in Mignone et al. (2009). 

Our simulation test results are very similar to the results using an HLL approximate Riemann 
solver scheme performed by Mignone et al. (2009) (see the lower panels of Fig. 12 in Mignone et 
al. 2009) except for their now smoother, when compared to the HLLD result, octagonal thin shell 
structure. We note that other previous relativistic and non-relativistic results for the MHD 2D 
rotor test problem using the HLL scheme (Del Zanna et al. 2003; Balsara & Spicer 1999) found a 
thin elliptical shell structure similar to our result in the xy plane (see the top-left panel in Fig. 5 of 
Del Zanna et al. 2003). Therefore the thin elliptical shell structure appears commonly developed in 
the rotor test problem if the HLL scheme is employed. Here we must speculate that the difference 
between the thin elliptical shell structure found by us and others using the HLL scheme, and 
the smoothed octagonal shell structure found by Mignone et al. (2009) using the HLL scheme is 
associated with a difference in the numerical diffusivity resulting from details of the piecewise linear 
reconstruction and constrained transport schemes used by Mignone et al. (2009). 

Figure 18 shows one-dimensional density profiles along the y and z axes for the different 
resolutions at time t = 0.4. In profiles along the y axis, the density in the central region is larger at 
lower resolution. Lower resolution gives a lower shell density and underestimates the propagation 
speed. Our simulation test results are very similar to those found using the HLL approximate 
Riemann solver scheme in Mignone et al. (2009), see the plus and star symbols of Fig. 14 in 
Mignone et al. (2009), except for the position of the thin shell. This difference is caused by the 
difference between the octagonal shape of the thin shell found by Mignone et al. and the elliptical 
shape of the thin shell found here. A protrusion associated with the octagonal-like thin shell along 
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the y axis extends beyond the smoother eUiptical shell. 

In summary, we obtain very similar results to previous studies of the 3D and 2D rotor test 
problem using the HLL approximate Ricmann solver scheme in RMHD. Wc conclude that our 
simulation code passes this test problem as successfully as other second-order codes. 
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Table 1. Models and Parameters 



Case 


a 


Vj/C 


Rj/a 


Pitch 


CPsa/2 


1.0 


0.2 


0.5 


constant 


CPsa 


1.0 


0.2 


1.0 


constant 


CPs2a 


1.0 


0.2 


2.0 


constant 


CPs4a 


1.0 


0.2 


4.0 


constant 


CPfa/2 


1.0 


0.3 


0.5 


constant 


CPfa 


1.0 


0.3 


1.0 


constant 


CPf2a 


1.0 


0.3 


2.0 


constant 


CPf4a 


1.0 


0.3 


4.0 


constant 


CPO 


1.0 


0.0 


0.0 


constant 


DPsa/2 


2.0 


0.2 


0.5 


decrease 


DPsa 


2.0 


0.2 


1.0 


decrease 


DPs2a 


2.0 


0.2 


2.0 


decrease 


DPs4a 


2.0 


0.2 


4.0 


decrease 


DPO 


2.0 


0.0 


0.0 


decrease 
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Table 2. Models, e-folding times and kink speeds 



Case 


Vj/c 


Rj/a 


Te/tc^ 




CPO 


0.0 


0.0 


3.75 


0.0 


CPsa/2 


0.2 


0.5 


3.75 


~ 


CPsa 


0.2 


1.0 


4.10 


0.06 


CPs2a 


0.2 


2.0 


4.35 


0.15 - 0.10^= 


CPs4a 


0.2 


4.0 


4.05 


0.20 - 0.16^^ 


DPO 


0.0 


0.0 


2.45 


0.0 


DPsa/2 


0.2 


0.5 


2.40 


~ 


DPsa 


0.2 


1.0 


2.45 


0.04 


DPs2a 


0.2 


2.0 


2.75 


0.12 - 0.06^ 


DPs4a 


0.2 


4.0 


2.70 


0.20 - 0.10^ 


CPO 


0.0 


0.0 


3.75 


0.0 


CPfa/2 


0.3 


0.5 


3.70 


< 0.05 


CPfa 


0.3 


1.0 


4.30 


0.12 - 0.06^ 


CPf2a 


0.3 


2.0 


5.60 


0.23 - 0.17^ 


cpria 


0.3 


-1.0 


5.30 


0.30 



^values indicated to nearest 0.05 and de- 
termined to ± 0.1 

'^values determined to it 5% 

'^values at early - late simulation times 



